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A field theoretic variational approach is introduced to study ion penetration into water-filled cylin- 
drical nanopores in equilibrium with a bulk reservoir (S. Buyukdagli, M. Manghi, and J. Palmeri, 
Phys. Rev. Lett., 105, 158103 (2010)). It is shown that an ion located in a neutral pore undergoes 
two opposing mechanisms: (i) a deformation of its surrounding ionic cloud of opposite charge, with 
respect to the reservoir, which increases the surface tension and tends to exclude ions form the pore, 
and (ii) an attractive contribution to the ion self-energy due to the increased screening with ion 
penetration of the repulsive image forces associated with the dielectric jump between the solvent 
and the pore wall. For pore radii around 1 nm and bulk concentrations lower than 0.2 mol/L, this 
mechanism leads to a first-order phase transition, similar to capillary "evaporation" , from an ionic- 
penetration state to an ionic-exclusion state. The discontinuous phase transition exists within the 
biological concentration range (~ 0.15 mol/L) for small enough membrane dielectric permittivities 
(em < 5). In the case of a weakly charged pore, counterion penetration exhibits a non-monotonic 
behavior and is characterized by two regimes: at low reservoir concentration or small pore radii, 
coions are excluded and counterions enter the pore due enforce electroneutrality; dielectric repulsion 
(image forces) remain strong and the counterion partition coefficient decreases with increasing reser- 
voir concentration up to a characteristic value; for larger reservoir concentrations, image forces are 
screened and the partition coefficient of counterions increases with the reservoir electrolyte concen- 
tration, as in the neutral pore case. Large surface charge densities (> 2 x 10~^ e/mr?) suppress the 
discontinuous transition by reducing the energy barrier for ion penetration and shifting the critical 
point towards very small pore sizes and molar reservoir concentrations. Our variational method is 
also compared to a previous self-consistent approach and yields important quantitative corrections. 
The role of the curvature of dielectric interfaces is highlighted by comparing ionic penetration into 
slit and cylindrical pores. Finally, a charge regulation model is introduced in order to explain the 
key effect of pH on ionic exclusion and explain the origin of observed time-dependent nanopore 
electric conductivity fiuctuations and their correlation with those of the pore surface charge. 

PACS numbers: 03.50.De,87.16.D-,68.15.-|-e 



I. INTRODUCTION 

Electrostatic forces induced by macroscopic dielectric 
bodies immersed in water regulate important phenom- 
ena such as the stability of colloidal suspensions [T], 
membrane assemblies , and ion selectivity by synthetic 
membranes |5 as well as in biological nanopores [4 . The 
image forces induced by the dielectric permittivity jump 
between dielectric bodies and the solvent surrounding 
them play a central role in these phenomena. Indeed, 
the equilibrium of similarly charged objects in a solvent 
is driven by the competition between the repulsive elec- 
trostatic interaction of their charges and the attractive 
van der Waals forces that originate from their low dielec- 
tric permittivity. This picture, valid at low ionic concen- 
trations, is also the basis of the DLVO theory [S]. A very 
similar competition is known to determine the perme- 
ability of biological channels that regulate ion exchange 
between the exterior and interior of cells. The strong 
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dielectric discontinuity between the membrane (dielec- 
tric permittivity ~ 2) and the water-filled channel 
(ew ~ 78) is at the origin of the high potential barrier for 
ion penetration into the pore. In 1969, Parsegian found 
that the energetic cost to move an ion from the bulk 
reservoir into a cylindrical pore of infinite length and ra- 
dius a = 0.2 nm is approximately 16 fcsT [1]. This result 
clearly suggests that at room temperature, any biological 
pore would be totally impermeable to ions. From the nu- 
merical solution of the Debye-Hiickel (DH) equation, it 
was later shown that the consideration of the finite length 
of the channel reduces this barrier up to 6 fc^T [7] for 
short enough pores. Levin [S] recently proposed an ap- 
proximative, but reasonably accurate, analytical solution 
of the DH equation for a finite length cylinder. He also 
showed that the introduction of surface charges can fur- 
ther reduce this energy barrier by attracting counterions 
electrostatically. Several perturbative approaches around 
weak-coupling (WC) theory [mean-field (MF) theory cor- 
rected up to one-loop or Debye-Huckel (DH) order] or 
strong coupling (SC) theory (virial expansion) [SI [TO] 
have been developed for charged slit pores [TTHT5] . con- 
centric cylinders [16] as well as for two like-charged di- 
electric cyhndrical bodies [17]. Because these methods 
tend to misevaluate ionic correlations, especially in the 
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presence of strong dielectric discontinuities, it becomes 
imperative to develop analytical or numerical tools valid 
over a larger parameter range. 

In the case of ions confined into a closed geometry char- 
acterized by a dielectric discontinuity, the existence of 
an infinite number of image charges significantly compli- 
cates reliable numerical calculations such as MC simu- 
lations [IH]- Although new MC algorithms for such sys- 
tems have been developed and applied to planar dielectric 
slabs [T9H23] . at the present, we are not aware of any ap- 
plication of these methods to dielectric cylinders where 
the curvature of the interface is likely to cause a further 
complication of the numerical task. When MC data for 
ionic penetration into dielectric cylinders become avail- 
able in the future, it will be interesting to test the va- 
lidity of approximative theoretical methods, such as the 
one presented in this article. 

Differential equations for electrostatic potentials, that 
partially take into account ionic correlations at a non- 
linear level (i.e. beyond the DH theory) were derived 
by Netz and Orland [21] within a field-theoretic varia- 
tional formalism. However these variational equations 
are too complicated to be solved analytically or even nu- 
merically [21]. A restricted tractable variational method 
for neutral single and double planar interface as well as 
for spherical systems was proposed in Ref. [25H27j . Be- 
cause of the piecewise inverse variational screening length 
introduced in 26J, the numerical minimization proce- 
dure is technically involved. We have recently proposed 
a simpler variational method to investigate the case of 
charged planar interfaces and slit pores [55] • The ap- 
proach is based on a generalized Onsager-Samaras ap- 
proximation |29j which assumes an electrostatic kernel 
with a uniform variational inverse screening length k„ 
that may differ from the bulk value. In this article, we 
extend this variational method to the case of neutral and 
charged cylindrical pores '30], by using a constant vari- 
ational Donnan potential 0o which enforces electroneu- 
trality in the charged pore [2S]. We show here that the 
extremization of the variational grand potential with re- 
spect to Ky and 00 yields two coupled variational equa- 
tions, 

/ dSa, = ypV(?,PM(e"*'^'^""^): (2) 

where we define the pore average as (^(r)) = ^^(r), 
Vp (p for pore) stands for the volume occupied by the con- 
fined ions, 5'p is the charged surface, Os is the uniform 
surface charge density (expressed in units of the elemen- 
tary charge e), qi denotes the ion valency, and p\,^i is the 
reservoir density of each ionic species. We also define the 
potential of mean force (PMF) 

$,(r;K„) = -ln^ = ^u;(r;K„) + g,0o (3) 
Ph,% 2 



which is the change in the excess electro-chemical poten- 
tial when an ion is brought from the bulk inside the pore 
at position r (in the following, all energies are rescaled by 
the thermal energy /cbT — l3~^). The quantity w{r, Ky) 
incorporates dielectric and electrostatic solvation forces 
(see Fig. [T]) and thus depends in a complex manner on 
Ky, and will be calculated in Section II. It is the Ky de- 
pendence of w that couples (t>o to without which (f)Q 
is simply the Donnan potential. The variational scheme 
employed here, which handles ionic correlations and im- 
age forces at a non-linear level, was recently applied to 
charged slit pores [55] • It was shown that the method i) is 
able to explore the regime between the WC and SC lim- 
its for charged single and double interfaces and ii) goes 
beyond the self-consistent methods used in nanofiltra- 
tion theories [351 131] • The ability of the present method 
to interpolate between the limiting WC and SC theo- 
ries makes it a valuable tool for exploring the physics of 
charged biological systems (lipid membranes, water chan- 
nels. . . ) where image forces play an extremely important 
role. Such a non-perturbative method is also essential 
for studying the possibility of an ionic liquid-vapor (L-V) 
transition in a nanopore, our principal goal here |30) . 
In order to clarify the need for such a method, we note 
that the for the weakly charged nanopores investigated 
here the V-phase can to a very good approximation be 
assimilated to a counter-ion only phase in the SC limit, 
because the despite the low values of the SC parameter, 
S ~ 0.1 the inequality (Eq. 4 of Jho et al, PRL 2008) 
delimiting the SC range of validity in a slab can, when 
extrapolated to a cylindrical nanopores of radius ~ 1 nm, 
still be satisfied thanks to the large value of the Gouy- 
Chapman length ~ 10 nm. Since the ionic L-phase is 
closer to the (salt) weak coupling limit, it is clear that a 
method is needed that spans the region between the two 
limiting laws. 

The variational equations obtained by Netz and Or- 
land [21] are equivalent to the closure equations estab- 
lished in the context of nanofiltration theories (see [32] 
I33j and |31j for a review). The closure equations were 
solved for spherical pores within a self-consistent approx- 
imation by Dresner [33], who observed a discontinuous 
phase transition from a high to a low ionic penetration 
state for decreasing reservoir concentration. Although 
the spherical geometry adopted by Dresner significantly 
simplifies the technical difficulties, it remains a toy-model 
since in reality, ion-penetration is controlled by thin chan- 
nels connecting the pore to the reservoir. The closure 
equations were later solved within the mid-point approx- 
imation by Yaroshchuk ^T] for a cylindrical pore. The 
mid-point approximation is equivalent to replacing the 
potential of mean force in the exponential of Eqs. 
by its value in the middle of the pore, which leads to an 
underestimation of repulsive image forces and solvation 
deficit effects. As will be shown below, this approxima- 
tion can overestimate the partition coefficients of ions in 
the pore by up to a factor of three. 

This paper is organized as follows. We introduce in 
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FIG. 1: Geometry for a cylindrical pore of infinite length and 
radius a. The cylindrical coordinates (r, 6, z) are defined in 
the figure. The ionic cloud is distorted due to the dielectric 
repulsion by the polarization surface charge density and the 
solvation deficit outside the pore. 



Section II the general lines of the field-theoretic approach 
for charged dielectric bodies and derive the variational 
grand potential for a cylindrical pore. Section III deals 
with ion penetration into a neutral pore of radius a. 
A first-order ionic exclusion transition from an ionic- 
penetration state to an ionic-exclusion state is found for 
a specific range of parameter values. We draw the phase 
diagram in the (a, pb) space. We then compare our re- 
sults with the self-consistent approach within the mid- 
point approximation [31] . In order to compare the ionic 
penetration into slit [5S] and cylindrical pores, we also 
compute the electrostatic potential and the variational 
free energy for two neutral concentric cylinders and study 
the limit of large cylinder radii with fixed separation. In 
Section IV, we first investigate the effect of a non-zero 
fixed surface charge on ion partitioning and show that 
the transition survives for sufficiently weak surface charge 
density, but disappears above a critical value because of 
the additional cost in electrostatic energy of the pore sur- 
face charge, leading to an enhanced ion penetration into 
the nanopore favorable to its screening. Finally, we de- 
velop our approach for the case where the surface charge 
density is regulated by the pH (surface charge regulation 
mechanism) which occurs experimentally for nanotubes 
in, e.g., PET membranes |34j where weak acid groups 
are located at the pore surface. The calculation of the 
Green's function for a cylindrical dielectric interface |35j 
is presented in Appendix A and for two concentric cylin- 
ders in Appendix B. 



II. MODEL 

In this section, we compute the variational grand po- 
tential of a symmetric electrolyte of dielectric permittiv- 
ity (identical to the bulk value) confined in a cylinder 
of radius a and length L. We consider in this article the 
limit of a sufficiently long cylinder (L ^ a) and neglect 
size effects. The cylindrical pore traversing a membrane 
is depicted in Fig.[T] The membrane of dielectric permit- 
tivity Em is salt free (i.e. k = 0) and the electrolyte is 



in contact with an external particle reservoir at the end 
boundaries of the cylinder. Hence the fugacity of ions is 
\i — e^*/Af, where A^ — h/y^2^7TmJ^BT is the thermal 
De Broglie wavelength of an ion i and /ij its chemical 
potential, and it is fixed inside the pore according to the 
chemical equilibrium condition, = Xi_b- 

The grand-canonical partition function of p species of 
interacting charges is 



p CO 



N,l 



2=nE 



The electrostatic interaction in Eq. Q is given by 
Ec = \ [ drdr'p(r)«c(r,r')p(r') 



(4) 



(5) 



where the charge distribution, in units of the elementary 
charge e, is 



P Ni 

i=i j=i 



(6) 



and /3s(r) ~ (Js6{r — a) is the negative fixed charge 
distribution (cts < 0) at the surface of the cylinder (1 
e nm~^ = 0.16 C m~^). The electrostatic potential 
Wc(r,r') is solution of 



where 



^;-i(r,r') = -^V[e(r)V<5(r-r')] 



e(r) = em6(r - a) + e„e(a - r) 



(7) 



(8) 



is the dielectric permittivity (where Q{x) is the Heavi- 
side distribution). Furthermore, the bulk self-energy of 
mobile ions that we substract from the total electrostatic 
energy m Eq. ^ IS 



(9) 



The bare Coulomb potential in the bulk is v^{r) = 
is/f, where the Bjerrum length is defined as £b — 
eV(47rewfcBr) ~ 0.7 nm in water at T = 300 K. After 
performing a Hubbard-Stratonovitch transformation and 
summing over Ni in Eq. Q, the grand-canonical parti- 
tion function becomes Q ~ J'Dcj) e~^^'^^/Zc, where the 
field-theoretic Hamiltonian is 



H\(b] = / dr 



[V(/.(r)]2-za(r)(/.(r)-^A,e' 



ao) 

where we have introduced the rescaled fugacities, — 

2 

A.e^^'cW. The average electrostatic potential -0(1") is re- 
lated to the fluctuating one 4>{r) by ■0(r) = i{(j){r)). The 
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factor Zc in the partition function subtracts the non- 
screened van der Waals contribution. Since we are exclu- 
sively interested in the salt-dependent part of the grand 
potential, we keep in the functional integral. The 
variational method consists in extremizing, with respect 
to the variational parameters, the first-order cumulant 
f2i, = f2o + ~ -f^o)o where the expectation values (. . .)o 
are evaluated with the variational Gaussian Hamiltonian 



[0(r) - *0o(r)] v^,\v, r') [0(r') - 



(11) 

and i7o — — 5trln(tio/fc)- The variational parameters 
are the Green's function i;o(r,r') and the electrostatic 
potential (/>o(r)- In the following, we consider the re- 



stricted case of a constant variational "Donnan" potential 
and uo(r,r') a solution of the inhomogeneous varia- 
tional Debye-Hiickel equation 

[-V(e(r)V) + e(r)K2(r)] wo(r, r') = ^e^b(v ~ r') (12) 

with a variational inverse screening length defined by 



k(v) = Ky Q{a — r). 



(13) 



The restricted variational choice made here is both sim- 
ple enough to lead to a tractable method and judi- 
cious enough to capture the essential physics for tight 
nanopores (a < 2 nm). After evaluating the functional 
integrals, the variational grand potential becomes 



= - ^ A, / exp 



'if 



{k^Ib - 5vq{v, r; k„)) - qi(t>n 



^ + '^f^ J <i£,(^dvo{r,r;Kyy^)~Svo{r,r;Ky)'^ + -as(l)o 



(14) 



where Vp — t:o?L and 5va{v,r\Ky) (see Fig. [2]) is the 
correction due to the presence of the nanopore to the 
variational Green function: 



-K„|r-r'| 



Vo(Y,Y ; Ky 



r — r' 



— -|-5wo(r,r';K„) (15) 



evaluated at r' = r [see Eq. (A9|], and is thus defined as 



AkY^ F,r,{k-Ky)ll,{K\v\) (16) 



m>0 



where we note = k'^ -\- and the prime on the sum- 
mation sign means that the term m = is multiplied 
by 1/2. The function Fm is a combination of modified 



Bessel functions /,„ and Km and is given in Eq. (AlOl. 
Note that we have Em for a — > oo and thus Svq 0. 

To find the physical meaning of the various contribu- 
tions in Eq. ( 14 ) , it is interesting to rewrite it in two 



ways. The classical thermodynamic equality in presence 
of surface effects, fly = —pVp + 'jSp with Sp = 2iTaL, 
allows us to separate the volumic bulk contribution, the 
pressure, which is independent of a (or Sp) 



P[Kv ) = 2^ Ai exp 



247r 



(17) 



and a surface contribution term, 7 ~ {fly + pVp) / Sp, 



^ ^ X^e'if 1^-^0/2 A _ g-gf5i;o(r,r;K„)/2-9i0o\ j_ "^t' 

2 ^ * \ / 16tt£b Jo 



d^(^6vQ{r,r;Kyy^) - 6vo{r,r; Hy)'^ + o-^'/'o (18) 



which is a function of nanopore characteristics (cts , a and 
ern via Svq), and vanishes for a — > 00. By maximizing 



the variational pressure in Eq. (17), we find the inverse 



screening length in the bulk, k^j given by the following 
implicit equation as a function of the fugacities A^: 



kI = AniB ^ QiK exp 



(19) 



Note that Eq. ( 19 ) leads to instabilities for large values of 



Xi if hard-core repulsion is not included Moreover, 



Eq. ( 19 1 has no solution for large fugacities such that 
qfKb(X)£B > 4. For low fugacities, Eq. ( [l9| yields the 
Debye-Hiickel Limiting Law (DHLL). Indeed, the bulk 
density is computed through 



Pj,b(A) = A, 



dp 



^ A, eT"""-^^ 



which allows to rewrite Eq. ([T9|) as 



Kl^4^reB}Q■p^,biX). 



(20) 



(21) 
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FIG. 2: Variation of Svo given in Eq. (161 with k„ for three 
values of r = 0,0.75, and 0.9 a (solid lines from bottom to 
top) and Em = 2, e„ = 78, and a = 0.84 nm. The dotted line 
corresponds to (f>o for |crs| = 5.4 X 10"** nm~'^ and pb — 0.1812 
mol/L. The left and right reference lines mark the solutions 
of just below and above the transition point, respectively. 



In the following, we will keep Ph i n the equations by re- 
placing Ae 2 '^>>'^b \yy jj^ (14). Note that to get p^. 



we calculate it from A, which is fixed in the grand canon- 
ical ensemble, using Eqs. (19)-(p0|. Solving Eq. (20) for 
/ii leads to the DHLL canonical ensemble result: 



ln(p,,bA3) - ^Kt,lB. 



(22) 



The second way to interpret Eq. (14 1 is by remark- 



ing that the first term on the rhs looks like the pres- 
sure of an ideal gas of ions "dressed" by their sur- 
rounding cloud, in a self-consistent or "external" field, 

2 

^(— K^j^s + <5wo(r, r; K„)) + Qi^o- It favors high Ky. The 
second and third terms are the correlation contribution, 
or the charging energy needed to create dressed ions [37] , 
modified by the nanopore. They come from "attractive" 
correlations between ions and their surrounding ionic 
cloud of opposite charge |37| . The second term is minus 
the DH pressure of a hypothetical bulk of inverse screen- 
ing length d the third one is a surface contribution 
due to the presence of the nanopore (dielectric exclusion 
and modified solvation effects). These two terms favor 
low ACi, because the third term, Q.d (Eq. A. 15), like the 
second bulk one, increases with k„, even if 6vQ{Ky) is 
a decreasing function of Ky (see Fig [2]). If, on the one 
hand, is set arbitrarily to (no screening, which cor- 
responds to "phantom ions" ) , one ends up with a simple 

barometric law, i.e. an ideal gas in the external potential, 

2 

^^uo(r, r; Q)+qi(j)o^ induced by the dielectric discontinu- 
ity at the pore surface. Within our restricted variational 
method this limit is the analog of the SC approximation 
(the ion-nanopore interactions dominate over the ion-ion 



correlations). If, on the other hand, we arbitrarily set 
Kv = Kb cosh[(/()o] > Kb and neglect the influence of the 
dielectric discontinuity on the Donnan potential (/jq, we 
obtain the analog of the WC PB/DH theory. As will be 
shown below by an application of our variational method, 
these simplifying limits are not sufficient for the task at 
hand. 

In summary, this system can be viewed as an ideal gas 
of dressed ions, whose accessible space is reduced by di- 
electric repulsion. Both the Boltzmann weight, which 
takes into account the feedback correlation associated 
with the fact that ions forming the surrounding cloud are 
themselves dressed ions, and the repulsive dielectric self- 
energy contribution depend on the variational screening 
parameter, Ky. The PMF, defined in Eq. ([s]), includes 
the potential w{r) that incorporates the solvation and 
image-charge interactions 



w{r) 



voir,r)-vliO) + Kb{X)£B 
{Kb - Ky)(.B + Svo{r, r; k^) 



(23) 



where Kb is defined in Eq. (19 1. The quantity q~w{r)/2 



is the difference between the excess chemical potential of 
ion i located at distance r in the nanopore and the excess 
chemical potential of the same ion in the bulk [28] . It is 
a decreasing function of Ky, since increasing k„ increases 
the electrostatic solvation gain of the hypothetic bulk 
(usual DH excess chemical potential term in —Kyis) and 
dvo{Ky) decreases with Ky (very slowly for very low and 
high Ky and abruptly at intermediate k.„ ~ l/(2a) , see 
Fig. [2]) , since the direct dielectric repulsion begins to be 
screened) . 



By minimizing Eq. ( 14 1 with respect to Ky, one exactly 
obtains Eq. ([T]) of the Introduction with 'w{r, Ky) defined 
in Eq. ( |23| . Note that for a bulk electrolyte (a — > 00) 
treated variationally, Eq. ([!]) yields the Debye-Hiickel re- 
lation Eq. (21) as long as the stability condition on fu- 
gacities \i,q~KbiB < 4 (pf, < 3 mol/L for monovalent 
ions) discussed in [^S] US' is satisfied. The drawback of 
working directly with Eq. ([l]) is that when it has three 
distinct solutions for Ky, one cannot distinguish between 
stable, metastable and unstable solutions. For this rea- 
son, in the following, we will analyze the pore-electrolyte 
model within the free energy minimization procedure and 
use Eq. ([I]) exclusively to show the importance of quan- 
titative errors induced by a previous self-consistent ap- 
proach [5T1 [55] . 

Finally, average ion concentrations in the pore are cal- 
culated using 



{P^{r)) 



-A. 



d{ny/Vp) 



(24) 

-Svo{r,r;K^)~qi4>o \ ^25) 



The ionic partition coefficients, which are measurable 
quantities, for instance in ion conductivity experiments. 
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are defined as 



(Mr)) 



drr-e-^"'('')-«"^°. (26) 



In terms of the A;^, the first (entropic) term in VL^, Eq. ( 14 1, 
becomes -J2zPi,bki- 

In the rest of the paper, we will consider a symmet- 
ric electrolyte, i.e. q+ = —q- — q. By differentiating 



the variational grand potential (14) with respect to 
one simply obtains the electroneutrality condition Eq. ([2]) 
which for a cylindrical nanopore is simply 



OS = qpbaTsinh{q<j)o) 
where we have defined the coefficient 



r = 



u{xa) 



(27) 



(28) 



which accounts for solvation and image corrections to 
mean- field theory (corresponding to F = 1). By inverting 
Eq. 1271) we find for cr^ < 



= - In 



qV/a 



-l + ^l + iqV/ay 



(29) 



where we have defined a — \as\/{pi,a) — Xjn/{2pi,), 
Xm = 2|(Ts|/a being the volume charge density of the 
pore. The average potential (f)Q increases in absolute 
value as 1/q \n\a/{2qT)] for a/qV S> 1. By injecting the 
solution Eq. ( 29 ) into the grand potential ( 14 ), we are left 



with a single variational parameter Ky that will be varied 
in order to find the optimal solutions to the variational 
problem. 



Q,/L 




(a) 
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(b) 

FIG. 3: Variational grand potential vs for tm — 2, e„ = 78 
and a = 0.84 nm. From bottom to top, bulk concentrations 
are (a) pt = 0.77 and 0.123 mol/L, and (b) pb = 0.181, 0.183, 
and 0.185 mol/L. 



III. NEUTRAL PORE 

In this section, we investigate the exclusion of ions 
from a neutral cylindrical pore. The electrolyte in the 
bulk reservoir is symmetric and composed of monova- 
lent ions {q = I). For a symmetric electrolyte the van- 
ishing surface charge, = 0, imposes 0o = through 
Eq. (27). In Fig. [3^ and [sjD are plotted the variational 
grand potential Qy vs the adimensional variational in- 
verse screening length Ibk^ for pore radius a = 0.84 nm, 
e„i — 2 (the case for lipid membranes), and various 
bulk concentration pi,. Figure [3^ shows that as one re- 
duces the reservoir concentration from pi, = 0.77 mol/L 
to pf, = 0.12 mol/L, the minimum of the grand poten- 
tial changes from Ibk^ ~ 1.1 ((p) = 0.2524 mol/L) to 
e^K^ = 0.037 {{p) = 3 X 10""* mol/L). In other words, 
the pore evolves from an ionic-penetration liquid state 
(L) to a quasi total ionic- exclusion vapor one (V). 

If one now slowly increases the bulk concentration from 
Pb = 0.12 mol/L to, for instance, pb = 0.181 mol/L, 
one notices the apparition of a second minimum at mod- 
erate IbKv — 0.34, shown in Fig. |3]d, which corre- 
sponds to a significant ion concentration in the pore. 



PL = 29 mmol/L. One sees that this minimum is 
metastable and the stable solution is the one correspond- 
ing to the ionic-exclusion state with pb = 0.44 mmol/L. 
If we keep increasing the reservoir concentration up to 
pI — 0.1832 mol/L, the values of both minima be- 
come equal, which indicates a phase coexistence: indeed, 
the equality of the two values of the grand potential, 
ily(Ky) = ^ly{Ky), indicates mechanical equilibrium be- 
tween both states. Finally, as pb is increased further 
to Pb = 0.185 mol/L, the ionic-exclusion state becomes 
metastable and the pore becomes penetrable to ions. 

This behaviour is the signature of a first- order phase 
transition over a certain parameter range for cylindrical 
nanopores. Note that such a transition does not appear 
to take place in slit-like pores although a continuous 
crossover from the presence to the absence of ions has 
already been observed for slit-like pores within a self- 
consistent calculation in Refs. [31] [33] and within the 
variational approach in Ref. [3S1[2H] (see Figs. 7-9). 

This discontinuous transition is characterized in Fig. [4| 
where we display for Cm = 2 the partition coefficients. 



given by Eq. ( 26 ), corresponding to the stable state of the 
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FIG. 4: Partition coefficients vs pb for t-m = 2, e^, — 78, and 
pore radii a — 1.96, 0.98 and 0.84 nm (solid lines, from top to 
bottom). The (red) circles display the metastable solutions, 
thus defining the coexistence window, red and red diamonds 
correspond to the solution of the self-consistent mid-point ap- 
proximation Eq. (31 1 for a = 1.96, and 0.98 nm. 



FIG. 5: Partition coefficient k = {p)/pb inside a neutral 
nanopore (q = 1, em — 2, = 78) vs the pore radius a for, 
from left to right, pt = 0.7,0.3,0.156,0.08 mol/L ("isobars"). 
Dotted (grey/red) lines show metastable (unstable) branches, 
light grey/red lines (guide for the eye) are the "boiling point" 
curve (bottom) and the "dew point" curve (top) and the dot 
is the critical point, pi = 0.075 mol/L, a* = 0.987 nm. 



variational grand potential, as a function of the reservoir 
concentration ph for three different radii. At large pore 
size a = 1.96 nni, the discontinuous jump is absent and 
the transition to the ionic-exclusion state is a continuous 
crossover. At a critical pore radius a* = 0.987 nm, a con- 
tinuous phase transition occurs, where a single minimum 
exists for Vl^ which evolves towards lower ion concen- 
trations in the pore with decreasing reservoir concentra- 
tion in a fast but continuous way (data not shown). For 
a — 0.98 nm (a < a*), the discontinuous jump at the co- 
existence bulk concentration exists, but with a very small 
jump for the partition coefficient. 

This is illustrated in Fig. [5] where is shown the par- 
tition coefhcient k versus the pore size for various pt,. 
The metastable and unstable branches are shown in red 
and the critical point, where the transition becomes sec- 
ond order, is for pi = 0.075 mol/L, a* = 0.987 nm. 
One notices the resemblance with the liquid-vapor (L- 
V) phase transition in bulk fluids, showing the "boiling 
point" curve and the "dew point" curve. Indeed, by in- 
terchanging the pore size, the partition coefficient of ions 
in the pore and the reservoir electrolyte density respec- 
tively with the temperature, the density and the pressure 
of a bulk liquid, one notices that the first-order ion ex- 
clusion transition in a nanopore becomes analogous to a 
bulk liquid- vapor (L-V) transition. 

To determine the parameter regime in which the dis- 
continuous transition takes place, the phase diagram is 
shown in Fig. [6] for — 78 and e„i = 1 to em = 4, 
where the lines correspond to the reservoir concentration 
at the coexistence point versus the pore radius, pl{a). 
Each curve corresponds to the coexistence line separat- 
ing the ionic-penetration state, L, (the area above the 
curve) from ionic-exclusion state, V (below the curve). 




0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 



a (nm) 

FIG. 6: Phase diagram characterizing the discontinuous phase 
transition for neutral pores {e^ = 78 and £2 = 1,2,3,4, from 
right to left). The critical line separates the ionic-penetration 
L state (area above the curve) from the ionic-exclusion V state 
(area below the curve) and ends at the critical point (dot) 
where the transition becomes continuous and then disappears. 



The lines end at the dependent critical point (a*,/?^), 
marked by a dot beyond which the transition disappears. 
One first notices that the parameter regime where the 
phase separation is observable is considerably reduced 
when increasing e^. The critical end-point evolves to- 
wards smaller pore sizes and higher reservoir densities. 
For ew = 78 and em from 1 to 4 (the curves in Fig. [6| , the 
critical pore size is a* — 1.44, 0.99, 0.77 and 0.65 nm while 
the critical reservoir concentration is pi = 17, 75, 221 
and 489 mmol/L, respectively. Higher values of em are 
not considered in this work since the coexistence line 
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would correspond in this case to very high bulk densities 
and extremely small pore sizes, where hard-core effects, 
which are not yet incorporated in our model, become 
non-negligible. 

It is clear from Fig.|6]that the membrane dielectric per- 
mittivity em, and thus dielectric repulsive forces, plays a 
central role in the mechanism of the phase transition. 
Moreover, it is known that the geometry influences the 
intensity of this dielectric repulsion. In order to inves- 
tigate curvature effects on the ionic exclusion mecha- 
nism, we compare in Fig. [7] the ion partition coefflcients 
in cylindrical (continuous lines) and slit pores (dashed 
lines). It is seen that as long as one is far from the tran- 
sition point, the ion density in a cylindrical pore of ra- 
dius a (diameter 2a) is very close to that in a slit pore 
of thickness d = a and the difference becomes smaller 
with increasing pore size a and/or reservoir concentra- 
tion pb. In order to explain this equivalence, we com- 
puted within the variational approach the partition coef- 
ficient of ions confined between neutral concentric cylin- 
ders of radius ai and 02 > ai (technical details are given 
in Appendix B). The result is illustrated in Fig. [t] for 
Pb = 0.075 mol/L and pb = 0.123 mol/L by two sets of 
points computed at a fixed distance between the cylin- 
ders, i.e. (5 = 02 — ai = 0.98 nm. The inner radius oi 
is varied from top to bottom between 0.7 and nm. At 
Pb = 0.123 mol/L and ai — 0.7 nm (the highest point), it 
is shown that the curvature of both cylinders is so weak 
that the ion concentration between the concentric cylin- 
ders is the same as that in the slit pore with d — S, a,s 
expected. By decreasing gradually ai from this value to 
0, one naturally recovers the ionic partition function in a 
cylindrical pore. For pb — 0.075 mol/L, the interpolation 
between the slit pore and the cylindrical pore occurs in 
a similar way, except that since the pore radius a — 6 
corresponds to the exclusion state, the partition coeffi- 
cient drops to nearly between ai = 0.119 nm (the sixth 
point from the top) and oi = 0.112 nm (the seventh point 
from the top). Hence, we conclude that the equivalence 
in ionic penetration for cylindrical pores of radius a and 
slit pores of thickness c? = a is simply the reminiscence of 
the concentric cylindrical case where for large values of 
the inner and outer radii (and the weak curvature of the 
dielectric interfaces) , the system qualitatively behaves as 
a slit pore. 

The existence of stronger dielectric forces in a cylindri- 
cal pore with respect to a slit one can thus be explained 
in terms of curvature. It is well-known that if the ion sees 
a dielectric interface with a positive curvature, i.e. if the 
interface is curved towards the ion, the experienced re- 
pulsion will be stronger than in the case where it would 
be placed face to a planar interface [38]. On the con- 
trary, if the ion is close to a negatively-curved interface, 
such as a protein or DNA, the dielectric exclusion will be 
weaker compared to a planar interface. This is the rea- 
son why the cavity correction factor of proteins measured 
in experiments is below one |39) . Therefore, the positive 
curvature of the pore is the essential ingredient for the 



Cyl, a = 1.96 nm 
Slit, d = a =1.96 nm 
Cyl, a = 0.98 nm 
Silt, d = a = 0.98 nm 
Cyl, a = 0.84 nm 

Con. cyl, a2-a1=0.9S nm, a1=0.0-0.7 nm 




Pb (mol/L) 



FIG. 7: Same figure as Fig. [4] for em ~ 1. Solid lines cor- 
respond to a cylindrical pore and dotted lines to a slit pore 
whose width is d = a. Solid circles correspond to a the parti- 
tion coefficient of an electrolyte confined between two concen- 
tric cylinders of a fixed width 5 — a2 — ai — a and show the 
crossover from a slit to cylindrical pore when ai varies from 
00 to 0. 



existence of a discontinuous phase transition. 



At the first-order variational level [see Eq. (14)], ion 
penetration is essentially driven by two opposing mech- 
anisms which contribute to the surface term 7, given in 
Eq. (18 1. On the one hand, due to "attractive" corre- 



lations between ions and their oppositely charged sur- 
rounding cloud, both the equivalent bulk term, K^/(247r), 
and the J d^ (Svoir, r, Kv^/S,) — Svo{r, r, Ky)) term 
vanish for Hy = and increase with Ky (the second term 
being sensitive to variations of Ky , because increasing Hy 
leads to the screening of Svq, see Fig. [2]). This cost in sur- 
face contribution is thus associated with the deformation 
of the ionic cloud of each ion due to dielectric repulsion 
and favors low Ky. 

On the other hand, for strong dielectric exclusion and 
weak surface charge, the presence of the pore walls de- 
creases the volume accessible to ions compared to an hy- 
pothetic bulk, and therefore the first (depletion) term in 
Eq. ( 18 ) contributes positively to the surface tension and 
decreases in magnitude with increasing Ky . When k„ 
increases, the dielectric repulsion becomes screened and 
the number of ions penetrating the nanopore increases 
(due to the Boltzmann law) , which decreases this surface 
cost. Indeed, the Gibbs adsorption equation at constant 
temperature is 



d7 = -{{p) - Pb)dn. 



(30) 



By increasing Ky, Eq. ( |23| ) leads to dp < 0, and since 
there is desorption in the nanopore ((p) < pb), dj < 0: 
this is favorable to the system and thus favors high Ky. 

Note that this liquid-vapor phase transition exists for 
bulk electrolytes but at very low temperature (T « 
100 K) for usual mineral salts BDH43I. The critical 
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temperature T* is thus shifted towards room temper- 
ature in pores of nanometer radius due to surface ef- 
fects. Indeed, essentiaUy due to dielectric repulsion, the 
low dielectric surface favors the vapor phase and thus 
shifts the first-order transition. We thus call this phe- 
nomenon capillary evaporation, just as for nano-capillary 
confined water, which undergoes condensation for hy- 
drophilic surfaces 4:4-47^ and evaporation for hydropho- 
bic ones gsmni. ^ 

The possibility for the existence of a first-order phase 
transition for electrolytes confined in spherical pores was 
first pointed-out by Dresner |33j . He used an approxi- 
mate Green's function, which is exact at the center of 
a spherical pore, and a self-consistent approach to show 
that the first-order phase transition survives as long as 
Cm «C Cw However, because of the numerous approxi- 
mations introduced in his work, Dresner |33j recognized 
that the quantitative predictions of this approach were 
not reliable. Yaroschuk [31^ later performed a similar 
self-consistent calculation using the correct Green's func- 
tion at the center of a cylindrical pore and showed that 
in this case a discontinuous phase transition to an ionic- 
exclusion state also coccurs. He also argued, however, 
that this phase transition is unphysical because it comes 
from the use of the DH equation, which ignores non- 
linear effects. By introducing a non-linear DH equa- 
tion (in analogy with the non- linear PB equation), he 
argued that the first-order phase transition should dis- 
appear. However, this PB-like generalization of the DH 
equation is incorrect as already argued by Onsager in the 
1930s [50J. Hence, we argue that, within the confined 
electrolyte model presented in this work, the existence of 
a first-order phase transition is physically sound. 

We finally analyze Yaroschuk's self-consistent (SC) ap- 
proach. It consists in finding self-consistently the screen- 
ing parameter k by defining as the average value of 
kI e"*^'"''^^ in the pore 31 . To pursue in a cylindrical 
geometry, he replaced the PMF within the integral by 
its value on the pore-axis, which in turn causes an un- 
derestimation of repulsive dielectric forces that become 
significantly stronger close to the pore wall. He obtained 
what we call the "mid-point approximation" for k: 

/t^^K^g-^Mo,.) (31) 

where Hb is defined in Eq. ( [2l| ) . The partition coefficients 
in a neutral pore obtained from the numerical solution 
of Eq. (31 1 are compared in Fig. [Ijwith the prediction of 
our variational approach for a — 1.96 and 0.98 nm. It is 
clearly seen that the mid-point approximation overesti- 
mates ionic penetration into the pore, which is due to the 
underestimation of dielectric forces, as stressed above. 
This point was also noticed for slit-pores in our previous 
work 28 . As one sees in Fig.|4j the error induced by the 
mid-point approach increases with decreasing pore size, 
which originates from the amplification of image forces 
when one decreases the pore radius. One notices that for 
a = 0.98 nm the partition coefficients as well as the po- 
sition of the transition point predicted by the mid-point 



approximation can deviate from the predictions of our 
variational method by 200-300 %. 



IV. CHARGED PORE 

A. Fixed surface charge 

In this part, we investigate ion penetration into a cylin- 
drical nanopore of non-zero negative fixed surface charge 
density, a-g < 0. We first note that the electroneutrality 



condition Eq. (27) can be written, using Eq. (261, as 



qpbO 



qpb 



(32) 



We introduce the good coion exclusion (GCE) limit which 
corresponds to the counterion-only case in the nanopore. 
In the case of slit pores, it was shown that the charge re- 
pulsion GCE limit is reached for small pore sizes, strong 
surface charge, and low bulk concentrations even with- 
out any dielectric effects [28] (F w 1, 2\as\/{qpbCi) ^ 1) 
and associated with the electric repulsion of colons. In 
the case of cylindrical nanopores, a second dielectric re- 
pulsion GCE limit exists for low \aa\ as soon as F ^ 
2\as\/{qpba). The partition coefficients of ions in the 
cylindrical pore become 



2H 
qpba 



and fc_ 



qaPb 
2|ctJ 



F^ < k. 



(33) 



The first equality shows that in the GCE regime, the av- 
erage counterion concentration in the pore, = k+pb, 
is independent of pb and depends only on the pore radius 
a and the surface charge a^. Hence, we are in a regime 
where dielectric repulsion plays only an indirect role (by 
excluding colons) in determining the average nanopore 
counterion concentration, which is solely determined by 
the global electroneutrality in the pore. 

We plot in Fig.[8]the partition coefficient of ions k± and 
the variational Donnan potential (po against the reservoir 
concentration pb for = 2, a = 0.84 nm and a weak sur- 
face charge, |crs| = 5.4 x 10"''' nm~^. One observes that 
at high bulk concentrations, k+ and fc_ decrease when pb 
decreases, as in the neutral case [and indeed (/)o ~ since 
qV ^ a in Eq. (29]), the slight difference fc+ — fc_ being 
given by Eq. ( 32 ijdown to a characteristic value pjj where 
a discontinuous jump towards a weak ionic-penetration 
state takes place. It is interesting to note that the Don- 
nan potential also exhibits a jump at this point. A 
first conclusion is the existence of a discontinuous tran- 
sition for charged pores at low \os\, with the coexistence 
value Pb((Ts) ^ Pbi'^s = 0). If one further decreases the 
reservoir concentration, fc_ evolves towards vanishingly 
small values while k+ abruptly changes by rapidly in- 
creasing. We thus reach the dielectric repulsion GCE 
regime, Eq. ( [SS] ), denoted by the dotted curve. Although 
the concentration of ions is low, it is not zero (in order to 
fulfill electroneutrality) and the vapor phase is thus now 
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FIG. 8: Partition coefficients of coions and counterions versus 
tlie reservoir concentration (above) and tlie effective Donnan 
potential (below) for em = 2, ~ 78, a = 0.84 nm and 



5.4 X 10"* nm" 



The dotted line in the top plot 



denotes the GCE regime Eq. ( 33 1 



a weak ionic-penetration phase. In this regime, \(f)o\ in- 
creases abruptly (as l/q\n[a/{2qr)]) because the dielec- 
tric exclusion dominates and thus gF <C a in Eq. ( 29 1 . 



Figure [9^ shows the partition coefHcients for a slightly 
stronger surface charge density \as\ = 1.35 x 1G~^ nm^^. 
In this case, the discontinuous jump disappears and the 
transition becomes continuous. Hence a large enough 
fixed surface charge density destroys the first-order phase 
transition. Figure|9]D shows the partition coefficients for a 
higher surface charge \as \ = 2 x 10~^ nm~^. In this case, 
not only the discontinuous phase transition, but also the 
turning point characterizing the abrupt change in the 
behaviour of disappears from the range of displayed 
Pb- In other words, above a characteristic surface charge 
and at low concentrations, the dielectric repulsion GCE 
regime sets in and the counterion penetration into the 
pore is solely determined by the global electroneutrality 
Eq. (|33l). 



In the previous section, we emphasized that the first- 
order nature of the transition is associated with the varia- 
tions of Svo with Ky, which influences both the nanopore 
modified DH correlations (which favor the low density 
V phase) and the entropic term (which favors the high 
density L phase). The absence of a first-order phase tran- 
sition when increasing as beyond a characteristic value 
is due to the increase of the surface contribution asso- 
ciated with the positive electrostatic energy of the pore 



surface charge (the last term in Eqs. (14 18)), which 



favors the higher ion concentration in the nanopore nec- 
essary for screening out the pore surface charge. The 
complete characterization of the phase transition is il- 
lustrated in Fig. [To] where is shown the phase diagram 
for several values of ct^. The main effect of the surface 
charge is to reduce the coexistence line and to shift the 
critical point towards smaller pore sizes and higher reser- 
voir densities. Comparison of Figs. 10 and [6] clearly shows 
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FIG. 9: Partition coefficients of coions and counterions versus 
the reservoir concentration for em = 2, tw = 78, a = 0.84 nm 
and \as\ = 1.35x10"^ nm"^ (above) and \as\ = 2x10"^ nm"^ 
(below). The dotted line in the top plot denotes the GCE 
regime Eq. ( 33 1. 



that the increase of the surface charge plays qualitatively 
the same role as an increase of the membrane dielectric 
permittivity. It is important to note that the smearing 
of the discontinuous phase transition takes place over a 
very narrow surface charge range: within the parameter 
regime considered in Fig. [lOj the transition completely 
disappears for \as\ > 6.0 x 10^'^ nm^^. 



It is well known that nature uses the very same pore 
surface chargemechanism to deal with strong image forces 
in water-filled ion channels [51^. As it was illustrated in 
the Section III, a neutral channel is impermeable to ions 
at low reservoir concentrations. In potassium channels, 
negatively charged carbonyl oxygens located at the pore 
surface [52l [53j reduce the potential barrier induced by 
image interactions and make the channel cation-selective. 
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FIG. 10: Phase diagram characterizing the discontinuous 
phase transition in charged pores for several values of the fixed 
surface charge. The critical line separates the ion-penetration 
state (area above the curve) from the ion-exclusion state (area 
below the curve) and ends at the critical point (dot), em = 2 
and Ew = 78. From top to bottom, the lines correspond to 
\as\ = 0,1,2,4,6 X 10"^ nm-^ 



B. Charge regulation mechanism 

In the previous part, we considered a boundary 
surface characterized by a fixed and uniform surface 
charge. However, in experiments on the conductivity of 
nanopores in PET membranes [34l[5j, the surface charge 
density is unknown and cannot be measured easily. If one 
wishes to make quantitative comparisons with these ex- 
periments, one should consider the fact that, in these ex- 
periments where the pore surface carries carboxylic acid 
groups, CTs increases with pH |55j . The pR in the reser- 
voir, defined as pH = — log^g pi, (H3O"''), is thus used as 
a control parameter. 

In this section, we introduce a charge regulation mech- 
anism where, in the presence of water, chemical groups 
located at the surface may become chemically active, 
leading to proton release from the acid groups and thus 
to a variable surface charge. As a first hint of the full 
problem, we use two approximations: i) we consider ex- 
clusively the case of trace hydronium ions, HaO"*", (spec- 
tator ions) with a bulk concentration (2 < pH < 12) 
significantly lower than that of salt ions; ii) the effective 
surface charge density, which becomes a function of the 
pH, is computed at the an effective mean-field level [56] . 
i.e. we consider that image forces, etc. for the hydro- 
nium ions in the pore and thermal fluctuations of the 
electric potential are implicitly accounted for via en ef- 
fective equilibrium constant. 

The charge dissociation mechanism at the surface is 
described by the following chemical equilibrium: 

[surface] - COOH + H2O ^ [surface] - COO" + H3O+ 

(34) 

characterized by the equilibrium constant Ka- Note the 



this effective Ka implicitly takes into account not only 
image forces but also all other non-electrostatic interac- 
tions at the pore surface. The source term [last term of 
Eq. (14)] is thus modified. In a simple Langmuir two- 
state model [Sni \57\ , each carboxylic group at the cylin- 
dric surface is a site, with surface density aQ, which can 
be empty, i.e. dissociated, with an energy per site —4>o, or 
occupied, i.e. associated with an hydronium ion, with an 
energy Ea which takes into account the chemical bound- 
ing and all other surface interactions and which, in a 
first approximation, is taken to be independent of k^. 
The grand potential for this adsorbed two-dimensional 
hydronium gas is then 



Q 



2D 



(35) 



where A+ = = Inp^^gO+e"'"'^''/^ is the fugacity of 
hydronium ions and Ea is the energy gained in the asso- 
ciated state. The equilibrium constant of Eq. (34) is thus 

_ TH3O+TCOO- _ p__E„+Kt 

7COOH 



defined as Ka 



* , where 7^ are 



activities coefficients of the relevant species. 
Instead of as 

M 



in Eq. (14), the source term becomes 



now 



2D 



S 



-(To In Q2D 



(36) 



= -O-O0O 



(7oln[l + 10f^^-P" e""^"] 



where pKa = —\ogiQ Ka- By minimizing the full grand 
potential with respect to (j)o , one obtains the electroneu- 



trality condition Eq. ( 27 ) with the effective surface charge 
density [58j 



A+ dn 



-(Jo 



2D 



S d\+ 



1 + lOPKa-pH g- 



(37) 



Figure [TT] illustrates for (To = 1 and 10 nm~^ the ionic 
partition coefficients (Fig. Ill) and the surface charge 
(Fig. [TTIj) versus the acidity of the solvent, pH — pKa. 
The low pH regime corresponds to the neutral pore limit, 
where hydronium ions of high bulk concentration almost 
totally neutralize the acid groups of the surface and the 
system is in the weak ionic penetration state. By in- 
creasing pH (or decreasing the hydronium concentration 
in the bulk) , the partition coefficient slowly increases un- 
til a characteristic value, pH'-^, beyond which one crosses 
the critical line and due to the high surface charge attrac- 
tion, the system evolves into the ionic penetration state. 
This transition is illustrated in Fig. [TT^ by the sudden 
jump of the coion and counterion partition coefficients, 
k±. The most remarkable prediction of our charge reg- 
ulation model is that the surface charge density displays 
a sharp discontinuous increase at the transition point. 
This behavior is obviously associated with the jump of 
the electrostatic potential (see Fig. [8| at the transition 
from the V- to the L-phase: salt ions induce the screening 
of the electrostatic potential in the pore 4>o ~ 0, which 
leads to a release of the spectator hydronium ions to the 
bulk and thus a sudden increase of the surface charge. 
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FIG. 11: (a) Partition coefficients of cations (biack lines) and 
anions (red fines) and (b) the reduced surface charge versus 
pH-pKa for em — 2, e^, = 78, pt ~ 1 mol/L and a — 0.617 nm. 
Sohd hues correspond to \ao\ = 1 nm~^ and dashed lines to 
(Jo — 10~^ nm~'^. The inset in the bottom plot illustrates the 
jump in the net surface charge for ao — 1 nm~^. 



The variations of k± with pH discussed above can 
explain the fluctuations found in experiments carried 
out on the conductivity of nanopores in PET mem- 
branes [311 I where a rapid switching between a high 
conductivity (HC) and low conductivity (LC) regime was 
observed. We argue that these fluctuations may corre- 
spond to the system being close to or at phase coexis- 
tence between both regimes, where the HC state would 
correspond to the liquid strong ionic penetration state 
and the LC state to the weak ionic penetration vapor 
state predicted within our nanopore model (we present a 
more detailed theoretical investigation of electric current 
fluctuations within our nanopore model elsewhere [30]). 
Furthermore, HC/LC switching is observed experimen- 
tally only within a narrow pH window and surface charge 
fluctuations are observed to be strongly correlated with 
conductivity ones [M|, both of which can be accounted 
for using our model, see Fig. [TT]b. and the inset of Fig. 
3b of [3(J] . Finally, it would be interesting to check exper- 



imentally our theoretical prediction for the ratio of the 
surface charge values, a^/aY — 0.51, at coexistence. 



V. CONCLUSION 

In this article, we developed a variational approach 
to investigate ion penetration into cylindrical nanopores. 
The variational scheme has already been developed in 
Ref. [21] and is based on a generalized Onsager-Samaras 
approximation which assumes a constant effective varia- 
tional screening length and a uniform Donnan potential 
enforcing elect roneutrality in the pore. These two physi- 
cally sound simplifications to the full variational problem 
allow us to handle the more complicated case of ions con- 
fined in a charged dielectric nanopore. 

In the first part of the work, we considered ionic ex- 
clusion from a neutral dielectric nanopore. It is shown 
that, when the electrolyte concentration in the reservoir 
or the pore size is decreased down to pb ~ 0.1 mol/L 
or a ~ 1 nm (see Fig. [6]), a first-order phase transition 
from an ionic-penetration state to an ionic-exclusion one 
emerges. The underlying physical mechanism this tran- 
sition is a competition between Gibbs adsorption which 
decreases the surface tension by favoring strong ion pene- 
tration and DH-type attractive correlations strengthened 
by the presence of the dielectric nanopore wall, which fa- 
vor low ion penetration. 

By comparing our predictions to the ones obtained us- 
ing a self-consistent mid-point approach f31|, it is shown 
that the latter can over-estimate ion concentration by a 
factor of 2 to 3. We also propose a continuous description 
from the slit geometry (thickness d), where no transition 
but rather a continuous crossover appears [55], to the 
cylindrical one (diameter 2c?), by investigating the case 
of concentric cylinders (separated by a distance d), thus 
highlighting the role of dielectric surface curvature. 

In the second part, we investigated the effect of a uni- 
form surface charge density on this ionic-exclusion phase 
transition. When varying bulk concentration, the coun- 
terion penetration into the cylindrical nanopore, induced 
by the charged surface in order to fulfill electroneutral- 
ity, exhibits a non-monotonic behaviour characterized by 
two regimes: i) a weak ionic-penetration (or dielectric 
repulsion GCE) regime, for low reservoir concentrations 
or small pore sizes, where the decrease of the colon par- 
tition coefficient with decreasing bulk salt concentration 
is even stronger than for a neutral pore, due to dielec- 
tric exclusion combined with electrostatic repulsion, and 
the counterion concentration remains low and is inde- 
pendent of pfj. ii) A second neutral-pore like ion pen- 
etration regime, reached above a characteristic reservoir 
concentration, where repulsive image forces are lower and 
the electrostatic potential is almost zero. The counter- 
ion partition coefficient fc+ increases with the reservoir 
concentration. For very small surface charge densities 
(< 2 X 10"'^ nm~^), the change between both regimes 
corresponds to a first-order phase transition. Increasing 
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the surface charge density shrinks the coexistence hne 
and moves the critical point towards smaher pore sizes 
and higher electrolyte concentrations. 

We finally consider the experimental case where the 
surface charge density varies with the bulk pH by intro- 
ducing a charge regulation mechanism. The first-order 
transition survives and qualitatively corresponds to ex- 
perimental results for conductivity fiuctuations measure- 
ments in nanopores [341 154) . Moreover, our charge regu- 
lation model, at an effective mean-field level, predicts a 
jump of the surface charge density as a funciton of pK at 
the transition, which is the signature of a sudden release 
of hydronium ions from the pore in going from the ionic 
Vapor to Liquid phase. 

Of course, our model has several limitations. First 
of all, due to the homogeneous form of the variational 
screening length, the dielectric/solvation deficit repul- 
sion is slightly underestimated. One interesting possible 
extension consists in introducing a piecewise screening 
length, as in Ref. but for an initial investigation 

of the type presented here the technical complications 
would outweigh the benefits. The model neglects also 
ion size, which becomes important for bulk systems at 
high concentration or strongly charged interfaces where 
ion concentrations may exceed the close packing value. 
For a cylindrical pore characterized by strong repulsive 
image forces where ionic concentrations are always far be- 
low their bulk value, we expect the ion size effect to play 
a less important role than in the bulk. However direct 
correlations, such as ion pairing, between "dressed" ions 
in the nanopore are neglected and could be taken into ac- 
count by extending the variational approach adopted in 
the present work to a second order cumulant expansion. 
Further investigation will be necessary to estimate the 
higher order contributions and we expect that it might 
change the parameter values at coexistence without mod- 
ifying our qualitative conclusions. 

Our present model also neglects ionic polarizability 
and treats the water solvent as a dielectric continuum 
of the same polarizability as in the bulk. It is important 
to note that due to various effects, the dielectric permit- 
tivity within the pore can be lower or higher than that of 
the bulk [SH]. The consideration of this additional com- 
plication requires a more detailed solvent model. We arc 
currently working on the generalization of our model to 
include the polarizability of charged ions, which should 
yield a more complete picture of the behaviour of large 
anions in confined geometries. The application of the 
present model to membrane nanofiltration with practical 
applications and direct comparison with experiments will 
be presented in a future work. 



Appendix A: Green's function in cylindrical 
coordinates : single cylinder 

In this appendix, we derive the Green's function 
vo{r,r') in the presence of a cylindrical dielectric in- 
terface of infinite length and radius a. The system 
is characterized by an electric discontinuity defined by 
e(r) — e<0(a — r) + e>0(r — a) and a Debye constant 
K{r) = K<0(a — r) + K>0(r — a), where r denotes the 
radial distance and ^ means r ^ a. We solve the Debye- 
Hiickel equation 

[-V(e(r)V) + e{r)K^{r)] vo{r, r') = I3e^6{r - r') (Al) 

by exploiting the cylindrical symmetry of the system, i.e. 
vo{r, r') = vo{z—z', 9—9', r, r'), where 9 stands for the az- 
imuthal angle. To this end, we expand wo(r, r') in Fourier 
space 



vo{r,r') 



+00 

E 



*m(e-8') 



+ CXD 



dk 

4^ 



^ik(z~z') 



VQ{r^ r'; m, k). 



(A2) 



By injecting this expansion into Eq. (Al) and using 
1 



S{9 



1 

2^ 



Sir - r')S{z - z')5{9 - 9') 

dk e^''-("-"') 



1 



+00 

27r ^ 

m— — oo 



(A3) 



Eq. ( Al I reduces to 



d'^vp ^ 1 dip 

Qj.2 J. Qj. 



4 + 7^ 



vp 



r 



(A4) 



= e 



where we have defined 

/3e^/(47re^). The homogeneous solutions of Eq. (A4) 



and = 



are the modified Bessel functions, vp{r,r']m,k) — 
A{r')I„i{>cr) + B{r')K„i{>cr). By taking into account 
the finiteness of the Green's function on the cylinder axis 
r = and at infinity r — > oo together with the continuity 
conditions 



(A5) 
(A6) 

(A7) 
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= ^p{t 
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one obtains the solution of Eq. ( Al ) in the form 
ir,r')=£^ + 5v< (r , r') . 



(A8) 
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The first term in ttie rhs is the usual DH potential in the bulk and the second term, which incorporates the dielectric 
discontinuity at the interface, reads for r' < r < a 



5v^{v,v') = — V cos[m(6' - 6*')] / dA;cos[/c(z - k<, K>)/„(>f<r)/„(>tf<r') 

TT — In 

m>0 



(A9) 



(AlO) 



and the prime on the summation sign means that the term m — must be multiplied by 1/2. Let us note that in the 
mid-point approximation (r = 0), the potential reduces to 



5v<{r = r' = 0,z = z',e = e') = ^ / 

7^ Jq 



e<>f<_ft'o(>f>a)^i(>f<a) — e>x>ifo(>«r<a)iiri(>ir>a) 



e<>f</<'o(>f>a)/i(>f<a) + e>x>/o(x:<a)i^i(>c>a) 
If the source is located outside the cylinder, one obtains in the region a < r < r' 

Sv^{r,r') = — cos[m{9 - e')] dkcos[k{z - z')]Gjn{k:, n^, Ky)Km{q>r)Krn{q>r') 

— In 

m>0 



Gm{k', K<, K>) — 



£>g>-^m(g<Q)C(g>Q) - ^<q<i-m{q>a)i'm{q<a) 

e<q<Km{q>a)I^(q<a) - eyqyI,n{q<a)K^{q>a)' 



(All) 

(A12) 
(A13) 



We also report the following coefficients introduced in Eq. (A15I 

Jo(fc;K) = /o(xra) - /i(xa) , 
Jmik;K) = /^(>ira) - /,„+i(xa)/„_i(xa) 



for TO > 0. 



(A14) 



Finally, in the calculation of the variational grand potential Eq. ( |14[ ), we perform analytically the integral over the 
radial distance r = |r| for the dielectric part which depends on Svq, which simplifies to 



m>0 



dk 



Fm{k;Kv\/£.,K> = 0)Jr,i{k;Kv\/£,) - Frn{k;Kv,Ky = 0)Jm{k;Kv) 



(A15) 



Appendix B: Green's function in cylindrical coordinates : concentric cylinders 



We report in this appendix the Green's function solution of Eq. ( A4 1 , for concentric cylinders of radius Oi and 02 
with oi < 02. We briefly explain in the end the computation of the variational free energy for this more complicated 
case. The system is characterized by a dielectric discontinuity defined by e(r) = e„i8(ai — r) + ewQ{r~' ai)0(a2 — r) + 
e?ri0('' — 0-2) and a Debye constant K(r) = KyQ{r — ai)Q(a2 — r) (we consider only the case of ions confined between 
the two cylinders). The solution of Eq. (|A4| satisfying the boundary conditions 
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(Bl) 
(B2) 

(B3) 
(B4) 
(B5) 



is given by 



g-K„|r-r'| 

i;o (r , r') = ^-b jr- + (5?;o (r, r') 



(B6) 
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with 



Svo{r,r') = —^y^ cos[m{9-e')] dkcos[k{z - z')] {Amik; fiy)I,ni^vr)Irn{>Cvr') 

TT — /n 

m>0 



(B7) 



where we have defined 



Bm{k', Ky) 
Cm{k] ^v) 



^m{k] Ky)hra{k] ^d) Cm{k] Ky)dm{k'^ ^d) 

dm{k] l^v)bm{k', Ky) 
^m{^': '^^))^m(^^ l^v) ^mik] ^v)dm{k] ^v) 

Cm{k; Ky)dm{k', Ky) 



^m{k] Ky)hyyi{k\ Ky) Cm{k] Ky)dm(k\ l^v) 

with the coefficients am{k]Ky), bm{k;Ky) and Cm{k; Ky) defined according to 

kKjn{xvai)I'^{kai) - ey,>iyIyn{kai)K'^{>Cyai) 
kIm{>Cya2)K'^{ka2) - e^„>f^X,„(fca2)/^(>fi,a2) 
Crn{k;Ky) = eyjKyK„iika2)K'^{>i:va2) - e„ikK„i{ka2)K'^{>i:ya2) 
dm{k]Ky) = eyjXyI„,{kai)I'jn{xyai) - trnklm{kai)l'^{xyai), 

and Ky = Ky + fc^. The Green's function evaluated at the same point is given by 

wo(r, r') = v^iO, 0) - Kylg + Svo{r, r') 



(B8) 



with 

(5wo(r,r) = — ^ / dfc {A„(fc; Ki,)/„(>f„r)^ + 5„(fc; K^^^i'ml^^i;?')^ + 2C„(fc; Kt,)/„(>f„r)fC™(>f„r)} . (B9) 



m>0 



The DH pressure from an hypothetical bulk solution is 
now given by 



^y 
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(BIO) 



and the spatial integrations in Eq. ( 14 1 are evaluated 



from r = fli to r = 02. As in the case with a single 



cylindrical pore, the obtained grand-potential energy fly 
is minimized with respect to Ky to find the optimal solu- 
tion to the variational problem. 
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